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Abstract 

Hydra  is  an  open-source,  platform-neutral  library  for  performing  Markov  Chain  Monte 
Carlo.  It  implements  the  logic  of  standard  MCMC  samplers  within  a  framework  designed  to 
be  easy  to  use,  extend,  and  integrate  with  other  software  tools.  In  this  paper,  we  describe 
the  problem  that  motivated  our  work,  outline  our  goals  for  the  Hydra  project,  and  describe 
the  current  features  of  the  Hydra  library.  We  then  provide  a  step-by-step  example  of 
using  Hydra  to  simulate  from  a  mixture  model  drawn  from  cancer  genetics,  first  using  a 
variable-at-a-time  Metropolis  sampler  and  then  a  Normal  Kernel  Coupler.  We  conclude  with 
a  discussion  of  future  directions  for  Hydra. 

Keywords:  Markov  Chain  Monte-Carlo,  Gibbs  Sampling,  Software  Library 
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1  Introduction 


Markov  Chain  Monte  Carlo  (MCMC)  is  a  method  of  performing  numerical  integration  on 
functions  that  can  expressed  as  distributions  (Metropolis  et  al.,  1953;  Hastings,  1970).  The 
strength  of  MCMC  is  that  it  can  simulate  from  distributions  without  requiring  the  densities 
to  be  properly  normalized.  This  makes  it  an  indispensible  tool  for  Bayesian  statistical  models, 
where  properly  normalizing  posterior  densities  is  often  impractical  or  impossible. 

After  an  initial  burn-in  period,  a  properly  constructed  MCMC  sampler  will  generate  a 
sequence  of  (non-independent)  samples,  Xo,Xi, . . .  ,Xj^  from  a  specihed  probability  distri¬ 
bution,  n.  Using  these  samples,  the  expectations  under  H  of  any  function  g  can  be  estimated 
by  computing  the  sample  path  average  E(g)  =  '^f=Qg{Xt). 

While  most  Markov  Chain  Monte  Carlo  algorithms  are  delightfully  simple,  there  are,  as 
of  this  writing,  only  two  software  packages  that  implement  general  MCMC  algorithms  for 
statistical  applications,  WinBUGS  and  FBM.  Both  of  these  have  important  limitations. 

WinBUGS  (Gilks  et  ah,  1992;  Gilks  et  ah,  1994b)  is  a  software  package  for  performing 
Bayesian  inference  using  Gibbs  sampling.  It  provides  tools  for  specifying  the  model,  running 
the  Gibbs  sampler,  and  monitoring  convergence  using  a  “point-and-click”  graphical  interface. 
A  noteworthy  feature  is  that  it  allows  models  to  be  specihed  using  either  a  text-based  notation 
or  a  graphical  model  created  with  the  DoodleBUGS  interface  (Spiegelhalter  et  al.,  1999). 

While  WinBUGS  is  mature  and  is  available  free  of  charge  from  the  MRC  Biostatistics 
Unit  web  site  (Stevens,  2000),  it  has  several  drawbacks.  First,  WinBUGS  is  designed  to  per¬ 
form  only  Gibbs  and  componentwise  Metropolis  sampling  and  does  not  allow  specihcation 
of  alternative  sampling  methods.  As  a  consequence,  WinBUGS  cannot  be  used  when  Gibbs 
sampling  or  Metropolis-within-Gibbs  are  inappropriate.  Second,  the  source  code  to  WinBUGS 
is  not  available  to  the  user.  This  makes  it  impossible  for  users  to  add  features  to  WinBUGS. 
Thus,  users  who  require  features  (such  as  statistical  distributions  or  sampling  methods)  not 
provided  by  WinBUGS  are  forced  to  abandon  the  package  entirely.  In  addition  the  inability  to 
access  the  source  code  prohibits  the  use  of  WinBUGS  for  experimentation  with  or  customiza¬ 
tion  of  sampling  algorithms.  This  prevents  WinBUGS  from  being  used  as  a  tool  for  research 
on  MCMC  methods. 

Radford  Neal’s  FBM  (“Flexible  Bayesian  Modeling”)  software  (Neal,  2000),  hrst  released 
in  1995,  is  a  less  well  known  package  that  implements  a  variety  of  MCMC  methods  and 
includes  the  C  source  code.  While  FBM  is  more  flexible  than  WinBUGS,  the  FBM  documentation 
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and  interface  are  considerably  more  difficult  to  understand. 

Both  WinBUGS  and  FBM  are  restricted  to  specific  operating  systems.  While  older  versions 
of  BUGS  were  available  for  Unix  systems,  the  current  version  is  available  only  for  systems 
running  versions  of  Microsoft  Windows.  FBM,  on  the  other  hand,  is  available  only  for  Unix 
systems.  In  addition,  neither  package  is  integrated  with  standard  statistical  tools.  This 
requires  the  user  to  learn  the  interface  of  an  additional  software  package  in  order  to  use 
MCMC. 

These  drawbacks  appear  to  have  discouraged  or  prevented  many  users  from  taking  advan¬ 
tage  of  the  considerable  effort  and  expertise  represented  by  WinBUGS  and  FBM.  As  evidence  of 
the  general  dissatisfaction  with  the  available  tools,  all  of  the  statisticians  we  have  observed 
using  or  researching  MCMC  write  their  own  custom  software.  This  results  in  considerable 
duplication  of  effort.  Worse,  since  properly  debugging  and  verifying  software  algorithms 
is  a  difficult  and  time-consuming  task,  it  is  likely  that  many  of  the  hand-written  software 
programs  contain  undetected  errors.  This  can  lead  to  the  presentation  of  faulty  analyses. 
Finally,  lack  of  integrated  software  support  for  MCMC  has  led  many  applied  researchers  to 
avoid  Bayesian  statistical  methods  entirely. 


2  Our  Approach 

Clearly,  there  is  a  need  for  better  MCMC  software.  Our  goal  is  is  to  produce  a  software  tool 
that 

1.  implements  standard  MCMC  techniques, 

2.  is  easy  to  use, 

3.  is  reliable, 

4.  is  applicable  to  a  wide  variety  of  problems, 

5.  allows  access  to  the  underlying  algorithms, 

6.  can  be  easily  customized  and  extended, 

7.  is  integrated  with  traditional  statistical  packages,  and 
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8.  is  available  on  all  common  computer  platforms. 

The  Hydra  MCMC  library  is  a  first  step  toward  providing  software  that  achieves  these 
goals.  Hydra  is  an  object-oriented  library  that  implements  the  logic  of  standard  MCMC 
methods.  Although  Hydra  can  be  used  directly  in  custom  MCMC  programs,  it  is  designed 
to  be  used  as  a  basis  for  MCMC  software  which  provides  a  more  user-friendly  interface  and 
which  is  integrated  with  standard  statistical  packages. 

Hydra  is  implemented  using  Java,  a  platform  independent  object-oriented  language  de¬ 
signed  for  general  programming  tasks.  We  selected  Java  (Joy  et  ah,  2000)  because  it  enabled 
the  library  to  meet  a  number  of  our  stated  goals.  First,  Java’s  support  of  formal  interfaces 
facilitated  the  construction  of  a  library  that  is  easy  to  use  without  sacrihcing  flexibility  and 
ease  of  extension.  In  particular,  the  use  of  interfaces  permits  users  to  reploace  exisiting  com¬ 
ponents  of  the  MCMC  algorithm  with  versions  which  are  tuned  to  the  specihc  problem.  This 
allows  the  user  to  extend  the  Hydra  package  to  support  new  problems  or  MCMC  techniques 
without  changing  the  existing  code.  Second,  Java  provides  features  that  reduce  common  pro¬ 
gramming  errors  and  is  supported  by  a  wealth  of  standard  libraries  and  programming  tools. 
Not  only  do  Java’s  features  make  it  easier  to  write  code  that  is  error-free,  but  they  also  make 
it  easier  to  locate  and  correct  bugs  that  do  exist.  This  supports  construction  of  a  reliable 
library  and  frees  time  otherwise  spent  debugging  for  the  development  of  additional  features. 
Third,  Java  provides  a  clear  interface  for  interacting  with  other  languages.  This  gives  a  well 
dehned  method  for  Hydra  to  be  used  with  existing  programming  languages  and  software 
applications.  Although  early  versions  of  Java  suffered  from  performance  problems,  recent 
versions  of  the  Java  virtual  machine  (runtime)  can  provide  speed  comparable  to  C  and  C++ 
for  numerical  computations  (Rijk,  2000;  Lewis,  2000;  Zachmann,  2000;  Schulman,  1997). 

The  remainder  of  this  text  assumes  a  basic  familiarity  with  the  Java  language  at  the  level 
of  Java  in  a  Nutshell  (Flanagan,  1997). 
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3  Constructing  Metropolis-Hastings  Samplers 


Using  Hydra 

Hydra  supports  the  full  generality  of  Markov  Chain  Monte  Carlo  by  providing  a  hierarchy  of 
classes  that  implement  the  most  common  MCMC  techniques.  Classes  implement  the  general 
Metropolis-Hastings  algorithm,  the  Metropolis  sampler,  and  the  Gibbs  sampler.  Hydra  also 
implements  the  multi-state  Adaptive  Metropolis  Sampling  (Gilks  &  Roberts,  1996)  algorithm 
that  forms  the  basis  of  Adaptive  Direction  Sampling  (Gilks  et  ah,  1994a)  and  Normal  Kernel 
Coupling  (Warnes,  2000).  We  will  focus  on  the  implementation  of  the  Metropolis-Hastings 
method  since  it  includes  the  others  as  special  cases. 

The  Metropolis-Hastings  algorithm  is  remarkably  simple.  Given  a  target  distribution 
of  interest  H,  corresponding  to  the  statistical  model,  an  initial  starting  location  Aq,  and  a 
proposal  distribution  Q{Xt),  each  iteration  of  the  sampler  consists  of  four  steps: 

1.  Propose  a  candidate  state  Y  using  the  proposal  distribution  Q{Xt),  which  may  depend 
on  the  current  state  Xf. 

Y  ^  Q{Xt) 

2.  Compute  the  Metropolis-Hastings  acceptance  probability 


a{Xt,Y) 


min 


min 


n(Y)  q(Y  ^  A',)  1 
’  ir(X,)  q(X,  ^  Y)  j 
p(Y)  q(Y  ^  X,)  ] 
’p(X,)  g{X,^Y)S 


where  tt  is  a  density  corresponding  to  the  target  distribution  H,  p{x)  oc  n{x)  is  an 
unnormalized  density,  and  q{Y  Xt)  =  q{Y\Xt)  is  the  conditional  density  of  Y  under 
Q{Xt). 
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3.  Accept  the  proposed  point  Y  and  set 


with  probability  a(Xt,Y),  otherwise 
Reject  the  proposed  point  and  set 


X, 


t+i 


X 


4.  Increment  time:  t  ^  t  +  1. 

The  sequence  of  X  values  generated  by  this  algorithm  converges  to  a  (dependent)  sample 
from  n  provided  the  proposal  distribution  Q  meets  certain  conditions  (Tierney,  1996). 

The  CustomMetropolisHastingsSampler  class  implements  the  logic  of  Metropolis-Hastings 
samplers  using  a  target  distribution  (model),  initial  state,  and  proposal  distribution  spec- 
ihed  by  the  user.  This  is  made  possible  by  requiring  the  objects  representing  the  target 
and  proposal  distributions  to  provide  certain  methods.  These  methods  are  dehned  by 
the  UnnormalizedDensity  and  GeneralProposal  interfaces,  respectively.  No  restriction 
is  placed  on  the  initial  state,  provided  it  is  compatible  with  the  user-specihed  target  and 
proposal  distributions. 

To  allow  flexible  reporting  of  the  progress  of  the  MCMC  sampler,  the  CustomMetropolis- 
HastingsSampler  maintains  a  list  of  user  dehned  objects  that  are  notihed  at  the  completion 
of  the  acceptance  step  of  each  iteration.  When  detailed  reporting  is  selected,  these  “listeners” 
receive  an  object  containing  a  great  deal  of  information  about  each  MCMC  iteration. 


3.1  The  UnnormalizedDensity  Interface  for  Target  Distributions 

Target  distributions  implement  the  UnnormalizedDensity  interface,  which  dehnes  two  meth¬ 
ods: 


public  double  unnormalizedPDF  (  Object  state  ); 
public  double  logUnnormalizedPDF  (  Object  state  ); 
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These  methods  compute  the  (log)  unnormalized  density  of  the  model  for  the  state  passed  as 
a  parameter. 


3.2  The  GeneralProposal  Interface  for  Proposal  Distributions 

Proposal  distributions  implement  the  GeneralProposal  interface,  which  has  4  methods: 


public  double  conditionalPDF  (  Object  next,  Object  current); 
public  double  logConditionalPDF  (  Object  next  ,  Object  current  )  ; 

public  double  transitionProbability  (  Object  from,  Object  to  ); 
public  double  logTransitionProbability  (  Object  from,  Object  to  ); 

The  methods  conditionalPDF  and  logConditionalPDF  compute  the  probability  of  gener¬ 
ating  the  object  next  when  the  current  state  is  current.  The  second  two  methods  perform 
the  same  computation,  but  reverse  the  arguments^. 


3.3  The  MCMCListener  Interface  for  Listener  Objects 

Objects  that  are  notified  at  the  completion  of  each  MCMC  iteration  implement  the  MCMCListener 
interface,  which  dehnes  one  method: 

public  void  notify  (  MCMCEvent  event  ); 

The  parameter  of  the  notify  method  is  an  object  containing  information  about  the  MCMC 
iteration.  This  information  can  be  used  by  the  object  in  various  ways.  Possibilities  include 
storing  the  current  state  to  a  hie,  displaying  it  on  a  plot,  and  computing  cumulative  statistics. 

When  detailed  reporting  is  disabled,  the  object  passed  to  notify  is  a  GenericChainStepEvent. 
This  object  has  a  single  held: 

public  Object  current  ; 

which  contains  the  current  state  (Xj)  of  the  sampler. 

When  detailed  reporting  is  enabled,  the  object  passed  to  notify  is  a  DetailChainStep- 
Event  which  has  the  additional  helds: 

^The  transitionProbability  and  logTransitionProbability  are  depreciated  and  will  not  be  required 
in  a  future  release  of  the  software. 


Table  1:  Intepretation  of  the  fields  of  the  DetailChainStepEvent 


Field 

Intepretation 

public 

Object  current; 

Xt 

public 

Object  proposed; 

F 

public 

double  proposedProb; 

piY) 

public 

Object  last; 

Xt-i 

public 

double  lastProb; 

p(Xi-i) 

public 

double  forwardProb; 

g(x|W_i) 

public 

double  reverseProb; 

q{Xt-i\Y) 

public 

double  probAccept; 

public 

double  acceptRand; 

uniform  value  used  to  accept /reject 

public 

boolean  accepted; 

was  Y  accepted? 

public 

double  acceptRate; 

average  value  of  Y) 

public 

Object 

proposed ; 

public 

Object 

last  ; 

public 

double 

lastProb  ; 

public 

double 

proposedProb  ; 

public 

double 

forwardProb  ; 

public 

double 

reverseProb  ; 

public 

double 

probAccept  ; 

public 

double 

acceptRand  ; 

public 

boolean 

accepted  ; 

public 

double 

acceptRate  ; 

These  helds  provide  a  great  deal  of  information  about  the  MCMC  iteration  and  are  useful 
for  debugging  and  for  evaluating  the  performance  of  different  proposal  distributions.  The 
interpretation  of  each  held  is  given  in  Table  1. 


4  Example 

The  classes  provided  by  Hydra  can  be  used  directly  in  compiled  Java  programs  or  interac¬ 
tively  with  various  Java-based  tools,  such  as  JPython  and  the  Omegahat  statistical  language. 
For  ease  of  presentation,  we  will  focus  on  the  pure  Java  interface. 

We  will  give  an  example  by  using  Hydra  to  construct  two  different  samplers  for  a 
Binomial-BetaBinomial  mixture  model  for  the  loss  of  genetic  material  in  esopageal  can- 
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Table  2:  Class  implementing  the  hierarchical  Binomial  Beta-Binomial  model  for  the  LOH 
data. 


package  org  .  omegahat  .  Simulation  .MIMC.  Examples  ; 
import  j  ava  .  lang  .  Math  ; 

import  org  .  omegahat .  GUtilities  .  ArrayTools  ; 

import  org  .  omegahat  .  Probability  .  Distributions  .  UnnormalizedDensity; 

public  class  Binomial_BetaBinomial_SimpleLikelihood  implements  UnnormalizedDensity 
{ 
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int  loh [ ] 

=  {  7, 

3, 

4,  3 

,  5  , 

4  , 

5  ,  3  , 

6  , 

1  2  , 

5  , 

3  ,  1  , 

3  , 

5  , 

3  ,  11,  2  , 

2  , 

2, 

3,  5, 

3, 

10 

4  , 

6  , 

3  ,  1  , 

4, 

5  , 

19,  5 

,  5, 

6  , 

5  , 

6  ,  2  , 

0  , 

0  , 

6,  4 

I; 

11 

int  n [ ] 

=  {17, 

15, 

17, 

18, 

15, 

,  15, 

15, 

19, 

16, 

15, 

18, 

19, 

18, 

19, 

19, 

2  1  , 

17, 

16  , 

12 

12, 

17, 

18, 

18, 

19, 

19, 

15, 

12, 

16, 

19, 

16, 

19, 

21  , 

15, 

13, 

20, 

16, 

17, 

13 

8, 

7, 

18, 

15}; 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 


/ /  unnormalized  binomial  density 

double  udb(int  x,  int  n,  double  pi)  {return  Math  .  pow  (  pi  ,  x )  *  Math  .  pow(  1 . 0— pi  ,  n— x);} 
/ /  unnormalized  beta— binomial  density 

double  udbb(int  x,  int  n,  double  pi  ,  double  omega)  ( 

int  r;  double  tmpC  =  1.0;  double  tmpl  =  1.0;  double  tmp2  =  1.0; 

for(r=0;  r  <=  (x  —  1  );  rH — h)  tmpO  *=  (  pi  +  ((double)  r)  *  omega); 

for(r=0;  r  <=  (n  —  x  —  1  );  rH — h)  tmpl  *=  (1.0  —  pi  +  ((double)  r)  *  omega); 

for(r=0;  r  <=  (n  —  1  );  rH — h)  tmp2  *=  (1.0  H-  ((double)  r)  *  omega); 

return  (  tmpO  *  tmpl  /  tmp2  )  ; 

} 

/ /  unnormalized  binomial— betabinomial  mixture  density 
double  ud_b_bb  (  int  x[],  int  n[],  double  eta, 

double  piO  ,  double  pil  ,  double  omegal )  ( 
double  retval  =  1.0; 

for(int  i=0;  i  <  x. length;  i+H-) 

retval  *=  (  eta)  *  udb  (x[i],  n[i],  piO)  H- 

(1.0  —  eta)  *  udbb  (x[i],  n[i],  pil,  omegal  )  ; 
return  retval ; 

} 

/ /  Constructor  / / 

public  Binomial_BetaBinomial_SimpleLikelihood  ()  (} 

/ /  Log  Unnormalized  Density  / / 

public  double  logUnnormalizedPDF  (  Object  parms )  ( 

return  Math  .  log  (  unnormalizedPDF  (  parms  ));} 

/ /  Unnormalized  Density  / / 

public  double  unnormalizedPDF  (  Object  paramObj  )  ( 
double  []  parms  =  ArrayTools  .  Otod  (  paramObj  ); 

double  et  a=parms  [  0  ]  ,  pi0=parms  [  1  ]  ,  pil=parms  [  2  ]  ,  omegal=parms  [  3  ]  ; 


/ /  check  range 

if  (  (eta<0.0)  ||  (pi0<0.0)  ||  (pil<0.0) 
(eta>1.0)  11  (pi0>1.0)  jj  (pil>1.0) 

return  0.0; 
else 

return  ud_b_bb(loh,  n,  eta,  piO  ,  pil 


(  omegal  <  0 . 0 )  | 
(omegal>0.5)  ) 


omegal  )  ; 
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cers.  We  first  show  how  to  implement  the  unnormalized  density  corresponding  to  Binomial- 
BetaBinomial  mixture  model  in  Java  so  that  it  can  be  used  with  Hydra.  Using  this  model 
we  create  a  variable-at-a-time  Metropolis  sampler,  and  then  a  Normal  Kernel  Coupler. 


4.1  Overview 

There  are  four  user-specihed  components  of  a  Metropolis-Hastings  sampler;  target  distri¬ 
bution  (model),  initial  state,  proposal  distribution,  and  uniform  random  number  generator. 
The  Hydra  library  provides  a  reliable  random  number  generator  and  a  selection  of  standard 
proposal  distributions,  leaving  the  user  to  construct  the  target  distribution  and  initial  state. 


4.2  Creating  a  Target  Distribution 

To  create  an  object  representing  the  target  distribution  (model),  the  user  needs  to  write  a 
Java  class  that  implements  the  UnnormalizedDensity  interface.  For  our  example  problem, 
we  wish  to  implement  a  class  for  the  Bayesian  hierarchical  model 

Xi  ~  ?7  Binomial (iVj,  TJ"!) 

-|-(1  —  77)  Beta-Binomial(iVj,  1^2) 

7]  ~  Unif(0, 1) 

TTi  ~  Unif(0, 1) 

7r2  ~  Unif(0, 1) 

UJ2  ~  Unif(0, 1/2) 


where  the  density  of  the  Beta-Binomial  distribution  is  dehned  as 


f{X,\N„n2,u;2) 


C) 


r(i;) _ ■~(.+g) 

r(^)r(i^)  r(n-a;+i^)r(n+^) 


A  Java  class  implementing  the  density  for  this  model  is  provided  in  table  2.  We  shall 
highlight  the  programming  details  that  allow  use  of  this  class  as  a  target  distribution. 

First,  we  need  to  indicate  to  the  Java  compiler  where  to  hnd  the  UnnormalizedDensity 
interface  that  this  class  will  implement.  This  is  accomplished  by  the  line: 

5  import  org  .  omegahat  .  Probability  .  Distributions  .UnnormalizedDensity; 
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Now  we  declare  the  class  and  indicate  that  it  implements  the  Unnormal izedDensity 
interface. 


7  public  class  Binomial_BetaBinomial_SimpleLikelihood  implements  UnnormalizedDensity 

Next,  the  class  needs  a  constrnctor  that  will  accomplish  any  reqnired  initialization,  such 
as  loading  the  observed  data.  In  this  case,  no  initialization  is  required  since  the  data  is 
hard-coded  into  the  class,  so  that  the  line 


40  public  Binomial_BetaBinomial_SimpleLikelihood  ()  {} 

is  sufficient. 

Now  our  class  must  provide  the  unnormal izedPDF  and  logUnnormalizedPDF  methods. 
These  methods  are  used  by  the  Metropolis-Hastings  sampler  to  compute  the  acceptance 
probability  for  a  proposed  state. 


42 

/ /  Log  Unnormalized  Density  / / 

43 

public  double  logUnnormalizedPDF  (  Object  parms ) 

{ 

44 

return  Math  .  log  (  unnormalizedPDF  (  parms  ));} 

45 

46 

/ /  Unnormalized  Density  / / 

47 

public  double  unnormalizedPDF  (  Object  paramObj  ) 

{ 

48 

double  []  parms  =  ArrayTools  .  Otod  (  paramObj  ); 

49 

double  et  a=parms  [  0  ]  ,  piO=parms  [  1  ]  ,  pil=parms[2 

]  ,  omegal=parms  [3]  ; 

50 

51 

//  check  range 

52 

if  (  (eta<0.0)  ||  (pi0<0.0)  ||  (pil<0.0)  ||  (omegal<0.0)  || 

53 

(eta>1.0)  11  (pi0>1.0)  jj  (pil>1.0)  jj  (omegal>0.5)  ) 

54 

return  0.0; 

55 

else 

56 

return  ud_b_bb(loh,  n,  eta,  piO  ,  pil  ,  omegal 

); 

57 

} 

In  this  example,  we  have  written  separate  functions  that  compute  the  unnormalized  density. 


so  these  methods  simply  convert  the  arguments  to  the  appropriate  type  (doubles),  check 
their  range,  call  the  appropriate  function. 

Note  that  the  interface  dehnes  the  argument  passed  to  the  unnormalizedPDF  and  logUnnormalized¬ 
PDF  as  an  Object.  The  user  must  decide  what  type  of  object  will  represent  the  model  param¬ 
eters.  For  most  purposes,  an  array  of  doubles  (double  [] )  is  an  appropriate  choice.  For  this 
reason,  the  predehned  proposal  methods  (see  Appendix  B)  all  operate  on  arrays  of  doubles. 

Any  other  type  of  object  may  be  used,  however,  this  will  require  the  user  to  implement  an 
appropriate  proposal  distribution. 
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4.3  Creating  a  Variable-at-a-time  Metropolis  Sampler 

Now  that  we  have  a  class  that  implements  the  unnormalized  density  for  the  Binomial- 
BetaBinomial  model,  we  can  construct  a  MCMC  sampler.  We  first  implement  a  variable- 
at-a-time  Metropolis  sampler.  The  complete  class  file  for  this  sampler  is  shown  in  table 

3. 

Table  3:  Class  implementing  a  variable-at-a-time  Metropolis  sampler  for  the  LOH  model. 


1  package  org  .  omegahat  .  Simulation  .MCMC.  Examples  ; 

2 

3  import  org  .  omegahat .  Simulation  .MCMC.  *  ; 

4  import  org  .  omegahat .  Simulation  .MCMC.  Proposals  ; 

5  import  org  .  omegahat .  Simulation  .MCMC.  Listeners  ; 

6  import  org  .  omegahat .  Simulation  .  RandomGenerators  ; 

7  import  org  .  omegahat .  Probability  .  Distributions  ; 

8 

9  public  class  Binomial_BetaBinomial_SimpleExample  { 

10  static  public  void  main  (  String  []  argv  )  throws  Throwable  { 

11 

12  CollingsPRNG  Administrator  a  =  new  GollingsPRNGAdministrator  ( )  ; 

13  PRNG  prng  =  new  CollingsPRNG  (  a  .  registerPRNGState  ( )  ); 

14 

15  UnnormalizedDensity  target  =  new  Binomial_BetaBinomial_SimpleLikelihood  ( ) ; 

16 

17  double[]  diagVar  =  new  double[]{  0.083,  0.083,  0.083,  0.042}; 

18 

19  SymmetricProposal  proposal  = 

20  new  NormalMetropolisComponentProposal  (  diagVar  ,  prng  ); 

21 

22  double[]  state  =  new  double  []{  0 . 9  0  ,  0.23,  0.71,  0.49  }; 

23 

24  CustomMetropolisHastingsSampler  mcmc  = 

25  new  CustomMetropolisHastingsSampler  (  state  ,  target,  proposal, 

26  pmg  9  true  )  ; 

27 

28  MCMCListener  1  =  new  ListenerPrinter  ( )  ; 

29  MCMCListenerHandle  Ih  =  mcmc.  registerListener  (  1  )  ; 

30 

31  mcmc.  iterate  (  10  ); 

32 

33  } 

34  } 


Again  we  provide  the  Java  compiler  with  the  locations  of  the  classes  we  will  be  using. 
This  time  there  are  hve  import  statements: 


import  org  .  omegahat  .  Simulation  .MCMC.  *  ; 
import  org  .  omegahat  .  Simulation  .MCMC.  Proposals.*; 
import  org  .  omegahat  .  Simulation  .MCMC.  Listeners  .*; 
import  org  .  omegahat  .  Simulation  .  RandomGenerators  .  *  ; 
import  org  .  omegahat  .  Probability  .  Distributions  .*  ; 
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After  declaring  the  object,  we  create  a  main  function  that  will  do  the  work  of  creating  and 
running  the  MCMC  sampler.  Within  main,  the  hrst  object  we  need  to  create  is  a  pseudo¬ 
random  number  generator.  The  Hydra  library  provides  an  implementation  of  the  Collings 
random  number  generator  (Collings,  1987),  which  can  be  created  using  the  2  lines: 


12  CollingsPRNG  Administrator  a  =  new  CollingsPRNGAdministrator  ( ) ; 

13  PRNG  prng  =  new  GollingsPRNG  (  a  .  registerPRNGState  ( )  ); 

Next,  we  need  to  instantiate  (create)  a  copy  of  our  class  that  implements  the  unnormalized 
density  of  the  model.  This  is  accomplished  by 


15  UnnormalizedDensity  target  =  new  Binomial_BetaBinomial_Likelihood  ( ) ; 

Now  we  instantiate  the  proposal  distribution.  For  a  Metropolis-Hastings  sampler,  there 
are  several  choices  (see  Appendix  B),  including  a  variable-at-a-time  random- walk  proposal 
using  a  normal  distribution.  This  is  implemented  by  the  NormalMetropolisComponentProposal 
class.  Its  constructor  allows  the  specihcation  of  a  proposal  variance  for  each  parameter.  We’ll 


use 

the  variance  of  the  parameters  under  the  prior: 

22 

double[]  diagVar  =  new  double[]{  0.083,  0.083,  0.083, 

0.042}; 

23 

24 

25 

SymmetricProposal  proposal  = 

new  NormalMetropolisGomponentProposal  (  diagVar  ,  prng 

); 

V  = 

Now  we  need  to  dehne  an  initial  state  for  the  sampler. 

:  0.90,  TTi  =  0.23,712  =  0.71,  a;2  =  0.49: 

We’ll  use  the  MLE,  which  is 

22 

double[]  state  =  new  double  []{  0 . 9  0  ,  0.23,  0.71,  0.49 

}; 

With  the  random  number  generator,  initial  state,  target  distribution,  and  the  proposal 
distribution  dehned,  we  can  create  the  actual  sampler: 

24 

25 

26 

GustomMetropolisHastings Sampler  mcmc  = 

new  GustomMetropolisHastingsSampler  (  state  ,  target, 

prng ,  true  )  ; 

proposal  , 

This  gives  us  a  working  MCMC  sampler.  The  hnal  parameter  is  an  optional  flag  indicating 
whether  the  sampler  should  report  all  of  the  details  about  the  MCMC  iteration  when  it  calls 
the  listeners,  or  whether  to  just  report  the  new  state.  We  wish  to  see  all  of  the  details,  so 
we  provide  the  value  true. 

We  need  to  attach  a  listener  to  the  MCMC  sampler  so  that  we  can  see  the  results  of  each 
iteration.  There  is  a  variety  of  predehned  listeners  (see  Appendix  C),  but  we’ll  start  with 
the  simplest  listener.  Its  notify  method  simply  displays  the  object  it  receives. 
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28  MCMCListener  1  =  new  ListenerPrinter  ( )  ; 

29  MCMCListenerHandle  Ih  =  mcmc.  registerListener  (  1  )  ; 

Finally,  with  the  MCMC  sampler  dehned  and  a  listener  attached,  we  are  ready  to  run 
the  MCMC  sampler.  This  is  accomplished  by  calling  the  MCMC  sampler’s  iterate  method 
with  the  number  of  iterations  to  perform: 


31  mcmc.  iterate  (  10  ); 


4.4  Running  the  Variable-at-a-time  Metropolis  Sampler 

When  combined  with  the  Hydra  library,  the  two  classes  we’ve  created  form  a  complete  Java 
program.  On  Unix-like  systems  with  the  standard  Sun  Java  tools  installed,  the  classes  can 
be  compiled  using  the  j  avac  command: 

>  javac  Binomial-BetaBinomial-SimpleLikelihood  .  java 

>  javac  Binomial_BetaBinomial_SimpleExample  .  j ava 

Once  the  classes  are  compiled,  the  MCMC  sampler  can  be  run  using  the  Java  interpreter 
by 

>  java  org  .  omegahat  .  Simulation  .MIMC.  Examples  .  Binomial_BetaBinomial_SimpleExample 

This  will  cause  the  MCMC  sampler  to  print  detailed  information  about  each  of  the  ten 
iterations  to  the  screen.  The  output  for  the  hrst  iteration  is: 


Chain  Step  Event  (with  details) 

Last  =  ContainerState  :  [  0.9  0.23  0.71  0.49  ] 

Last  Prob  =  -359.046964566765 

Proposed  State  =  ContainerState:  [  0.9  0.4751068415506061  0.71  0.49  ] 
Proposed  Prob  =  -423.26454869568283 

Current  State  =  ContainerState:  [  0.9  0.23  0.71  0.49  ] 

Forward  Prob  =  -0.0369493853787235 
Reverse  Prob  =  -0.0369493853787235 
Acceptance  Prob  =  -64.21758412891785 
Acceptance  Val  =  0.658405257229882 
Accepted?  =  false 

Acceptance  Rate  =  0.0 


This  gives  the  current  and  proposed  states,  the  value  of  the  unnormalized  density,  the  forward 
and  reverse  proposal  probabilities,  the  acceptance  probability,  a  flag  indicating  whether  or 
not  the  proposed  state  was  accepted,  the  new  state,  and  the  cumulative  acceptance  rate. 
Note  that  the  unnormalized  density  and  probabilities  are  reported  on  the  log  scale.  The 
output  for  all  10  iterations  is  given  in  Appendix  D 
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4.5  Enhancing  the  Variable-at-a-time  Metropolis  Sampler 

This  example  can  be  enhanced  in  a  number  of  ways.  First,  the  class  can  be  modihed  to  use  a 
different  proposal  method.  To  use  a  (complete-state)  random-walk  Metropolis  sampler,  sim¬ 
ply  replace  the  NormalMetropolisComponentProposal  with  a  NormalMetropolisProposal. 
Alternatively,  the  user  could  dehne  a  custom  proposal  distribution  and  use  it  instead. 

Second,  it  is  impractical  to  store  and  interpret  all  of  the  detailed  information  produced 
by  using  the  StepListenerPrinter  listener  for  more  than  a  few  iterations.  Instead,  we 
would  like  to  store  just  the  current  state  to  a  disk  hie.  This  is  accomplished  by  replacing 
the  StepListenerPrinter  object  with  a  StrippedListenerWriter.  Change  lines  28  and 
29  to 


The  1 .  close  0  ;  command  makes  sure  that  the  hie  that  is  used  to  store  the  MCMC  output 
is  properly  closed  once  the  MCMC  iterations  are  complete. 


Now  that  the  output  is  being  stored  to  a  disk  hie,  it  is  reasonable  to  increase  the  number 
of  iterations.  Naturally,  this  is  done  by  changing  the  value  in  the  mcmc.  iterate  call  to  the 
desired  value,  say  10,000. 

Compiling  and  running  the  modihed  class  now  generates  a  data  hie  containing  10,000 
MCMC  iterations.  This  output  can  be  read  into  a  standard  statistical  package  for  compu¬ 
tation  of  diagnostics  and  to  perform  inference.  For  this,  we  have  found  the  CODA^  package 
of  MCMC  diagnostics,  which  exists  in  versions  for  both  R  and  S-PLUS,  particularly  helpful. 
For  either  version,  the  commands 

library ( coda ) 

mcmc  .data  <—  mcmc(  as  .  matrix  (  read  .  table  ( ”MCM0.  output  ”  ) )  ; 

will  load  the  CODA  library  (provided  it  is  installed)  and  properly  import  the  MCMC  data. 
A  selection  of  diagnostics,  plots,  and  summaries  is  then  available.  For  instance,  the  default 
CODA  plots  and  summary  statistics  for  our  10,000  iterations  are  shown  in  Figure  1. 

^See  Appendix  A. 3  for  information  on  obtaining  CODA. 
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Iterations 


Figure  1:  CODA  plots  for 


Density  of  VI 

N  =  10000  Bandwidth  =  0.01267 

Density  of  V2 

N  =  1 0000  Bandwidth  =  0.003038 

Density  of  V3 

N  =  10000  Bandwidth  =  0.02709 

Density  of  V4 

).0  0.2 

0.4  0.6  0.8  1.( 

N  =  30000  Bandwidth  =  0.03909 

Table  4:  Class  implementing  a  Normal  Kernel  Coupler  for  the  LOH  model. 


1  package  org  .  omegahat  .  Simulation  .MIMC.  Examples  ; 

2 

3  import  org  .  omegahat .  Simulation  .MCMO.  *  ; 

4  import  org  .  omegahat .  Simulation  .MCMO.  Proposals  ; 

5  import  org  .  omegahat .  Simulation  .MCMO.  Listeners  ; 

6  import  org  .  omegahat .  Simulation  .  RandomGenerators  ; 

7  import  org  .  omegahat .  Probability  .  Distributions  ; 
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32 

33 
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35 
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45 

46 
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public  class  Binomial_BetaBinomial_SimpleExample_NKC  { 

static  public  void  main  (  String  []  argv  )  throws  Throwable  { 

CollingsPRNG  Administrator  a  =  new  GollingsPRNGAdministrator  ( ) ; 

PRNG  prng  =  new  CollingsPRNG  (  a .  registerPRNGState  ( )  ); 

UnnormalizedDensity  target  =  new  Binomial_BetaBinomial_SimpleLikelihood  ( ) ; 

double[][]  Var  =  {{  0.003,  0.0  ,  0.0  ,  0.0  }, 

{  0.0  ,  0.001,  0.0  ,  0.0  }, 

{  0.0  ,  0.0  ,  0.012,  0.0  }, 

{  0.0  ,  0.0  ,  0.0  ,  0.007}}; 

HastingsCoupledProposal  proposal  =  new  NormalKernelProposal  ( Var ,  prng  ); 
int  numComponents  =  200; 

MultiDoubleState  stateO  =  new  MultiDoubleState  (  numComponents  ); 
for  (int  i=0;  i  <  numComponents  /  2  ;  i++) 

state0.add(  new  double  []  {  0 . 9  0  3  ,  0.228,  0.708,  0.486  }  ); 

for  (int  i=numComponents  /  2  ;  i  <  numComponents;  i++) 

state0.add(  new  double  [  ]  {  0 . 0  7  8  ,  0.831,  0.230,  4.5e  — 9  }  ); 

Custom  Hastings  Co  up  led  Sampler  mcmc  = 

mcmc  =  new  CustomHastingsCoupledSampler  (  stateO  ,  numComponents, 

target ,  proposal ,  prng, 
false  )  ; 

ThinningProxyListener  pL  =  new  ThinningProxyList ener  ( numComponents  )  ; 
MCMCListenerHandle  pLh  =  mcmc.  registerListener  (pL  ) ; 

MCMCListenerWriter  1  =  new  StrippedListenerWriter  (”NKC.  output”  ) ; 
MCMCListenerHandle  Ih  =  pL  .  registerListener  ( 11  )  ; 

mcmc .iterate(  10000  ); 

1 . close  ( )  ; 

} 

} 
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4.6  Implementing  the  Normal  Kernel  Coupler 


To  implement  the  Normal  Kernel  Coupler  (NKC)  introduced  by  Warnes  (2000),  only  a  three 
changes  need  to  be  made  to  our  example  class.  First,  we  use  a  different  proposal  distribution. 
Second,  we  initialize  a  set  of  initial  values  rather  than  a  single  value.  Third,  we  use  the 
CustomHastingsCoupledSampler  class  instead  of  the  CustomMetropolisHastingsSampler 
class.  The  complete  source  code  for  the  modihed  class  is  given  in  table  4. 

The  proposal  distribution  for  the  NKC  is  implemented  by  the  class  NormalKernel- 
Proposal.  Its  constructor  requires  two  arguments,  a  random  number  generator,  and  a 
matrix  that  specihes  the  variance  for  the  normal  kernel.  For  our  example,  the  proposal  is 
instantiated  by  the  lines: 
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double  [  ]  [  ]  Var  =  {  { 

0. 

003, 

0.0  , 

0 

.0 

0.0 

18 

{ 

0. 

0  , 

0.001, 

0 

.0 

0.0 

}, 

19 

{ 

0. 

0  , 

0.0 

0. 

.012, 

0.0 

}, 

20 

{ 

0. 

0  , 

0.0 

0 

.0  , 

o 

o 

o 

7}}; 

21 

22 

HastingsCoupledPropO! 

3al 

pro] 

posal  = 

new  No 

rmalKernelProposal  ( Var  ,  prng  ); 

The  NKC  maintains  a  set  of  current  states  that  must  be  initialized.  We  use  a  MultiDoubleState, 
which  holds  a  list  of  double  values,  to  hold  the  initial  values. 


24  int  numComponents  =  200; 

25 

26  MultiDoubleState  stateO  =  new  MultiDoubleState  (  numComponents  ); 

27  for  (int  i=0;  i  <  numComponents  /  2 ;  i+H-) 

28  state0.add(  new  double  []  {  0 . 9  0  3  ,  0.228,  0.708,  0.486  }  ); 

29 

30  for  (int  i=numComponents  /  2 ;  i  <  numComponents;  i+H-) 

31  state0.add(  new  double  [  ]  {  0 . 0  7  8  ,  0.831,  0.230,  4.5e  — 9  }  ); 


In  this  case,  we’ve  initialized  half  of  the  values  to  each  of  the  two  local  maxima. 

The  logic  of  multi-state  MCMC  samplers  is  implemented  by  the  CustomHastingsCoupled- 
Sampler  class.  This  class  is  instantiated  using  6  parameters,  the  set  of  initial  states,  the 
number  of  current  states  to  maintain,  the  target  (model)  distribution,  the  proposal  distri¬ 
bution,  a  random  number  generator,  and  a  flag  indicating  whether  to  report  the  details  of 
the  iteration: 


33  CustomHastingsCoupledSampler  memo  = 

34  mcme  =  new  CustomHastingsCoupledSampler  (  stateO  ,  numComponents, 

35  target  ,  proposal  ,  prng  , 

36  true  )  ; 


Although  we  could  have  simply  used  a  StrippedListenerWriter  this  would  generate  a 


very  large  output  hie  by  writing  out  the  entire  set  of  200  current  states  at  each  iteration. 
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Instead,  we  use  a  ThinningProxyListener  class,  which  “thins”  the  events  it  receives  by  a 
specihed  factor  before  passing  them  on: 


38  ThinningProxyListener  pL  =  new  ThinningProxyListener  (numComponents )  ; 

39  MCMCListenerHandle  pLh  =  mcmc.  registerListener  (pL  ) ; 

40 

41  MCMCListenerWriter  1  =  new  StrippedListenerWriter  (”NKC.  output”  ) ; 

42  MCMCListenerHandle  Ih  =  pL  .  registerListener  ( 11  ) ; 


Running  the  sampler  now  will  output  the  complete  state  once  every  200  iterations. 


5  Conclusions  and  Future  Directions 

Our  example  has  shown  that  the  Hydra  MCMC  library  makes  it  easy  to  create  differ¬ 
ent  Metropolis-Hastings  samplers  without  extensive  programming.  This  should  encourage 
additional  statisticians  to  experiment  with  and  use  the  Metropolis-Hastings  method. 

We  hope  that  the  Hydra  library  will  form  the  basis  of  a  set  of  MCMC  tools  that 
are  easy  to  use,  robust,  and  complete.  In  particular  we  intend  to  integrate  Hydra  with 
the  statistical  tools  R,  Spins,  and  SAS,  as  well  as  the  new  Omegahat  statistical  computing 
language  (Temple  Lang,  2000;  Chambers,  2000;  Bates  et  ah,  2000).  These  interfaces  promise 
to  provide  flexible  and  powerful  interactive  environments  for  MCMC. 

Other  goals  for  the  Hydra  library  include 

•  visual  tools  for  specifying  and  monitoring  MCMC  simulations 

•  support  for  distributed/parallel  computing 

•  a  library  of  target  distributions  corresponding  to  common  statistical  models,  such  as 
GLM’s  and  mixture  models. 


20 


References 


Bates,  D.,  Chambers,  J.,  Cook,  D.,  Dalgaard,  P.,  Gentleman,  R.,  Hornik, 
K.,  Ihaka,  R.,  Leisch,  F.,  Lumley,  T.,  M  achler,  M.,  Masarotto,  G., 
Murrell,  P.,  Narasimhan,  B.,  Ripley,  B.,  Sawitzki,  G.,  Temple  Lang,  D., 
Tierney,  L.,  &  Venables,  B.  (2000).  The  Omega  project  for  statistical  computing, 
web  site,  http://www.omegahat.org/. 

Chambers,  J.  M.  (2000).  Users,  programmers,  and  statistical  software.  Journal  of 
Computational  and  Graphical  Statistics  . 

COLLINGS,  B.  J.  (1987).  Compound  random  number  generators.  Journal  of  the 
American  Statistical  Association  82,  525-527. 

Flanagan,  D.  (1997).  Java  in  a  Nutshell.  O’Reilly  &  Associates,  second  edition. 

ClLKS,  W.  R.  &  Roberts,  G.  O.  (1996).  Strategies  for  improving  MCMC.  In  Markov 
Chain  Monte  Carlo  in  Practice,  pages  89-114.  Chapman  &  Hall. 

Gilks,  W.  R.,  Roberts,  G.  O.,  &  George,  E.  I.  (1994a).  Adaptive  direction 
sampling.  The  Statistician  43,  179-189. 

Gilks,  W.  R.,  Thomas,  A.,  &  Spiegelhalter,  D.  J.  (1992).  Software  for  the  Gibbs 
sampler.  In  Computing  Science  and  Statistics.  Proceedings  of  the  2frd  Symposium  on  the 
Interface,  pages  439-448.  Interface  Foundation  of  North  America  (Fairfax  Station,  VA). 

Gilks,  W.  R.,  Thomas,  A.,  &  Spiegelhalter,  D.  J.  (1994b).  A  language  and 
program  for  complex  Bayesian  modeling.  The  Statistician  43,  169-177. 

Hastings,  W.  K.  (1970).  Monte  Garlo  sampling  methods  using  Markov  chains  and  their 
applications.  Biometrika  57,  97-109. 

Joy,  B.,  Steele,  G.,  Gosling,  J.,  &  Bragha,  G.  (2000).  The  Java  Language 
Specification.  Addison-Wesley,  second  edition,  also  at 
http : // j  ava . sun . com/ docs/books/ j Is. 

Lewis,  J.  P.  (2000).  Java  versus  G/G++  benchmarks,  web  site, 
http : / / WWW . idiom . com/~zilla/ Computer/ j  avaCbenchmark . html. 


21 


Metropolis,  N.,  Rosenbluth,  A.  W.,  Rosenbluth,  M.  N.,  &  Teller,  A.  H. 
(1953).  Equations  of  state  calculations  by  fast  computing  machine.  Journal  of  Chemical 
Physics  21,  1087-1091. 

Neal,  R.  (2000).  Software  for  flexible  Bayesian  modeling  and  Markov  chain  sampling. 
Web  Site,  http://www.cs.toronto.edu/~radford/fbni.software.htnil. 

Rijk,  C.  (2000).  Binaries  vs  byte-codes.  Ace’s  Hardware  (web  site) 

http : //www . aceshardware . com/Spades/read . php?article_id=153. 

ScHLFLMAN,  A.  (1997).  Java  on  the  fly.  Wehreview.com  (web  site) 
http : //webreview.  coni/wr/pub/97/07/25/grok/. 

Spiegelhalter,  D.  J.,  Thomas,  A.,  &  Best,  N.  G.  (1999).  WinBUGS  Version  1.2 
User  Manual.  MRC  Biostatistics  Unit. 

Stevens,  A.  (2000).  The  BUGS  project.  Web  Site, 
http : //www.mrc-bsu. cam. ac .uk/bugs/welcome . shtml. 

Temple  Lang,  D.  (2000).  The  Omega  project:  New  possibilities  for  statistical  software. 
Journal  of  Computational  and  Graphical  Statistics  . 

Tierney,  L.  (1996).  Introduction  to  general  state-space  markov  chain  theory.  In  Markov 
Chain  Monte  Carlo  in  Practice,  pages  59-74.  Chapman  &  Hall. 

Warnes,  G.  R.  (2000).  The  Normal  Kernel  Coupler:  An  adaptive  Markov  Chain  Monte 
Carlo  method  for  efficiently  sampling  from  multi-modal  distributions.  PhD  thesis. 
University  of  Washington. 

Zaghmann,  G.  (2000).  Java/C++  benchmark,  web  site, 
http : //www . igd . f hg . de/~zach/benchmarks/. 


A  Installing  Hydra 

The  Hydra  Java  package  is  available  in  two  forms,  as  a  Jar  hie  (Hydra,  jar)  containing 
only  the  compiled  classes  and  as  a  gzipped  tar  hie  (Hydra,  current  .tar  .gz)  containing  the 
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full  source  code  as  well  as  the  compiled  classes.  Both  hies  are  available  from  the  Hydra 
web  page  located  at  http://www.warnes.net/Hydra. 


A.l  Installing  the  jar  File 

Download  Hydra,  jar  and  append  the  full  path  to  the  jar  hie  to  the  CLASSPATH.  For 
example,  if  Hydra,  jar  has  been  placed  in  the  directory  /home/user/ jars/,  the  proper 
command  for  setting  the  CLASSPATH  using  sh  compatible  shells  is 

>  CLASSPATH=$CLASSPATH :  /  home  /  user/jars/Hydra.jar 

>  export  CLASSPATH 

and  using  csh  compatible  shells 

>  setenv  CLASSPATH  SCLASSPATH :/ home/ user  /  j  ar  s /Hydra  .  j  ar 


A. 2  Installing  the  Full  Source 

Download  Hydra,  current  .tar  .gz.  It  can  then  be  unpacked  using  GNU  tar  via 

>  tar  — xvzf  Hydra . current . tar . gz 

which  will  unpack  a  directory  tree  with  root  “Hydra”.  The  Java  hies  and  source  code  are 
contained  in  directories  under  Hydra/org/omegahat 

The  location  of  this  directory  then  needs  to  be  added  to  the  Java  class  path.  If  the  direc¬ 
tory  tree  was  unpacked  in  /home/user/jsrc/  this  can  be  accomplished  using  sh  compatible 
shells  by 

>  CLASSPATH=$CLASSPATH :  /  home  /user/jsrc/Hydra 

>  export  CLASSPATH 

and  using  csh  compatible  shells 

>  setenv  CLASSPATH  SCLASSPATH:  /  home/ user  /  j  src /Hydra 


A. 3  Other  Packages 

Two  additional  Java  packages  may  be  required  to  use  particular  features  of  the  Hydra 
library.  Visual  Numerics’  JNL  and  Omegahat.  In  addition,  the  CODA  package,  in  conjnnction 


23 


with  either  R  or  SPLUS  statistical  packages,  provides  a  useful  suite  of  tools  for  evaluating  and 
making  inference  using  MCMC  output. 


•  Visual  Numerics’  JNL  library  is  required  for  several  of  the  Hydra  classes,  in  particular 
those  used  in  the  Binomial-BetaBinomial  example  given  below.  JNL  is  available  free  of 
charge  from  the  Visual  Numerics  web  site: 

http : //www. vni . com/products/wpd/jnl/. 

•  The  Hydra  MCMC  classes  were  designed  to  be  compatible  with  the  Omegahat  statisti¬ 
cal  programming  system.  Omegahat  provides  an  interactive  environment  for  statistical 
programming  and  analysis  and  is  under  active  development  by  the  Omegahat  project, 
http : //www . omegahat . org. 

•  The  statistical  package  R  is  a  free  re-implementation  of  the  S  language  and  may  be 
obtained  free  of  charge  from  http :  //www .  r-proj  ect .  org. 

•  The  CODA  package  of  MCMC  diagnostics  and  other  tools  for  Splus  can  be  obtained 
from 

http : //www.mrc-bsu. cam.ac.uk/bugs/classic/coda04/readme . shtml. 

A  version  for  R  can  be  obtained  from  http://www-fis.iarc.fr/coda. 
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B  Predefined  Proposal  Distributions 


Hydra  provides  a  selection  of  predefined  proposal  methods  for  the  Metropolis,  Metropolis- 
Hastings,  and  Hastings-Coupled  techniques. 


Metropolis  Samplers 


Class  Name 

Description 

NormalMetropolisProposal 

NormalMetropolisComponentProposal 

normal  random-walk  proposal 
variable-at-a-time  random-walk 
proposal 

Metropolis-Hastings  Samplers 


Class  Name 

Description 

NormalProposal 

(hxed)  normal  proposal 

NormalMetropolisProposal 

normal  random-walk  proposal 

NormalMetropolisComponentProposal 

variable-at-a-time  random-walk 
proposal 

MixtureProposal 

finite  mixture  proposal  using 
specihed  components 
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Hastings-Coupled  (Multi-Chain)  Samplers 


Class  Name 

IndependentHastingsCoupledProposal 

Adapt iveNormalMetropolisProposal 

Adapt iveNormalProposal 

NormalKernelProposal 
Adapt iveNormalKernelProposal 

LocallyAdaptiveNormalKernelProposal 

KernelDirectionSampler 


Description 

Wrapper  for  independent 
Metropolis-Hastings  Samplers 
(Variance)  Adaptive  Normal 
Metropolis  Proposal 
(Mean,  Variance)  Adaptive 
Normal  Proposal 
Normal  Kernel  Coupler 
(Variance)  Adaptive  Normal 
Kernel  Coupler 

Locally-Adaptive  Normal  Kernel 
Coupler 

Kernel  Direction  Sampler 
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C  Predefined  Listeners 


A  variety  of  predefined  listeners  are  available.  These  allow  monitoring  various  features  of 
the  MCMC  simulation  and  give  several  storage  methods. 


Class  Name 

Description 

Accept  ance  Writer 

Stores  the  cumulative  acceptance  rate  to  a  file 

C  ovari  anceWriter 

Stores  the  cumulative  covariance  matrix  to  a  file 

DistanceListener 

Computes  the  observed  and  expected  acceptance  rate, 
step  distance,  step  distance  conditional  on  acceptance 

Distance  Writer 

Stores  the  observed  and  expected  acceptance  rate, 
step  distance,  step  distance  conditional  on  acceptance 
to  a  file 

Histogram  Writer 

Stores  a  cumulative  histogram  of  the  current  states  a  file 

ListenerGzipWriter 

Stores  the  current  state  to  GZIP  compressed  file 

ListenerPrinter 

Prints  the  event  passed  to  notify() 

ListenerWriter 

Stores  the  event  passed  to  notify()  to  a  file 

Mean  Writer 

Stores  the  cumulative  mean  vector  to  a  file 

PosteriorProb  Writer 

Stores  the  (unnormalized)  posterior  probability  of  the 
current  state  to  a  file 

Quantile  Writer 

Stores  the  cumulative  quantiles  to  a  file 

StepListenerPrinter 

Prints  MCMCStepEvents 

StrippedListenerGzip  Writer 

Stores  the  current  state  to  a  GZIP  compressed  file 

StrippedListener  Writer 

Stores  the  current  state  to  a  file 

ThinningProxyListener 

A  proxy  for  other  listeners  that  thins  the  reported  events 
by  a  specified  factor,  eg  1  out  of  every  100 
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D  Output  from  Binomial  BetaBinomial  Example,  java 


>  java  org  .  omegahat  .  Simulation  .IVOvIC!.  Examples  .  Binomial_BetaBinomial_Example 
Chain  Step  Event  (with  details) 

Last  =  ContainerState  :  [  0.9  0.23  0.71  0.49  ] 

Last  Prob  =  -359.046964566765 

Proposed  State  =  ContainerState:  [  0.9  0.4751068415506061  0.71  0.49  ] 

Proposed  Prob  =  -423.26454869568283 

Current  State  =  ContainerState:  [  0.9  0.23  0.71  0.49  ] 

Forward  Prob  =  -0.0369493853787235 

Reverse  Prob  =  -0.0369493853787235 

Acceptance  Prob  =  -64.21758412891785 
Acceptance  Val  =  0.658405257229882 

Accepted?  =  false 

Acceptance  Rate  =  0.0 

Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.71  0.49  ] 

Last  Prob  =  -359.046964566765 

Proposed  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.49  ] 

Proposed  Prob  =  -360.0165066537818 

Current  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.49  ] 

Forward  Prob  =  -0.06327351354448152 

Reverse  Prob  =  -0.06327351354448152 

Acceptance  Prob  =  -0.969542087016805 
Acceptance  Val  =  0.31056162077494043 
Accepted?  =  true 

Acceptance  Rate  =  0.5 

Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.49  ] 

Last  Prob  =  -360.0165066537818 

Proposed  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 
Proposed  Prob  =  -360.1951841070053 

Current  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 

Forward  Prob  =  0.6154440530408976 

Reverse  Prob  =  0.6154440530408976 

Acceptance  Prob  =  -0.17867745322348583 

Acceptance  Val  =  0.13272539253007873 

Accepted?  =  true 

Acceptance  Rate  =  0.6666666666666666 

Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 

Last  Prob  =  -360.1951841070053 

Proposed  State  =  ContainerState:  [  1.02451456  2  4463  2  77  0.23  0.45610096830883196 
0.4225189574152196  ] 

Proposed  Prob  =— Infinity 

Current  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 

Forward  Prob  =  0.23049155040119496 

Reverse  Prob  =  0.23049155040119496 

Acceptance  Prob  =  —Infinity 

Acceptance  Val  =  0.335025050367706 

Accepted?  =  false 

Acceptance  Rate  =  0.5 


Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 


Last  Prob  =  -360.1951841070053 

Proposed  State  =  ContainerState  :  [  0.9  0.1856011046441249  0.45610096830883196 
0.4225189574152196  ] 

Proposed  Prob  =  -363.1065248979661 

Current  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 

Forward  Prob  =  0.3116872397632934 

Reverse  Prob  =  0.3116872397632934 

Acceptance  Prob  =  -2.9113407909608213 

Acceptance  Val  =  0.14726731467399157 

Accepted?  =  false 

Acceptance  Rate  =  0.4 

Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 

Last  Prob  =  -360.1951841070053 

Proposed  State  =  ContainerState:  [  0.9  0.23  0.1825256980628881  0.4225189574152196  ] 
Proposed  Prob  =  -364.9924350935826 

Current  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 

Forward  Prob  =  -0.1255457772139429 

Reverse  Prob  =  -0.1255457772139429 

Acceptance  Prob  =  -4.797250986577353 

Acceptance  Val  =  0.4310649477090058 

Accepted?  =  false 

Acceptance  Rate  =  0.3333333333333333 

Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.4225189574152196  ] 

Last  Prob  =  -360.1951841070053 

Proposed  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.5879703533000055  ] 
Proposed  Prob  =  -359.85442835072524 

Current  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.5879703533000055  ] 

Forward  Prob  =  0.3415983954458079 

Reverse  Prob  =  0.3415983954458079 

Acceptance  Prob  =  0.0 

Acceptance  Val  =  0.36058319097411967 

Accepted?  =  true 

Acceptance  Rate  =  0.42857142857142855 

Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.5879703533000055  ] 

Last  Prob  =  -359.85442835072524 

Proposed  State  =  ContainerState:  [  0.6  69  65  8  70  69  5  04394  0.23  0.45610096830883196 

0.5879703533000055  ] 

Proposed  Prob  =  -361.2666064185844 

Current  State  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.5879703533000055  ] 

Forward  Prob  =  0.00517213125315924 

Reverse  Prob  =  0.00517213125315924 

Acceptance  Prob  =  -1.412178067859145 

Acceptance  Val  =  0.9268299028867995 

Accepted?  =  false 

Acceptance  Rate  =  0.375 

Chain  Step  Event  (with  details) 

Last  =  ContainerState:  [  0.9  0.23  0.45610096830883196  0.5879703533000055  ] 

Last  Prob  =  -359.85442835072524 

Proposed  State  =  ContainerState:  [  0.9  0.2099774552506214  0.45610096830883196 

0.5879703533000055  ] 

Proposed  Prob  =  -360.42346255905113 

Current  State  =  ContainerState:  [  0.9  0.2099774552506214  0.45610096830883196 
0.5879703533000055  ] 

Forward  Prob  =  0.32110939780366615 
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Reverse  Prob  =  0.32110939780366615 
Acceptance  Prob  =  -0.5690342083258884 
Acceptance  Val  =  0.3311132575995816 
Accepted?  =  true 

Acceptance  Rate  =  0.4444444444444444 

Chain  Step  Event  (with  details) 

Last  =  ContainerState :  [  0.9  0.2099774552506214  0.45610096830883196 

0.5879703533000055  ] 

Last  Prob  =  -360.42346255905113 

Proposed  State  =  ContainerState:  [  0.9  0.2099774552506214  0.15833697164158184 
0.5879703533000055  ] 

Proposed  Prob  =  -365.44521171685466 

Current  State  =  ContainerState:  [  0.9  0.2099774552506214  0.45610096830883196 
0.5879703533000055  ] 

Forward  Prob  =  -0.2084655958574132 
Reverse  Prob  =  -0.2084655958574132 
Acceptance  Prob  =  -5.0217491578035265 
Acceptance  Val  =  0.7418864833851747 
Accepted?  =  false 

Acceptance  Rate  =  0.4 
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